Este documento es la memoria de nuestra práctica para la asignatura Aprendizaje Computacional de la Mención en Computación del Grado de Informática de la Universidad de Murcia.
El proyecto consiste en el estudio de diferentes modelos de aprendizaje automático. Para ello, se ha utilizado la base de datos Credit Approval de UC Irvine Machine Learning Repository, que contiene información sobre la concesión o denegación de créditos bancarios.
En primer lugar, se ha llevado a cabo un análisis exploratorio de la base de datos. Hemos identificado atributos numéricos y categóricos y hemos renombrado las columnas según la información que hemos encontrado en un repositorio de Kaggle. Hemos proseguido con una serie de análisis monovariable y multivariable discriminando entre créditos concedidos y denegados. Esto nos ha permitido hacernos una idea de la relación entre los diferentes atributos. Para cerrar esta sección, hemos probado un Análisis de Componentes Principales (PCA) para tener una referencia sencilla con la que comparar los modelos de aprendizaje.
Después del análisis, hemos probado diversos modelos de clasificación supervisada utilizando el lenguaje R y la librería caret, aplicando técnicas de preprocesado, ajuste de hiperparámetros y evaluación cruzada. El dataset se ha dividido en datos de entrenamiento y datos de test para comprobar la eficacia de los modelos. En total, se han probado cuatro algoritmos representativos de distintos paradigmas de aprendizaje automático.
Comenzamos instalando las librerías de R necesarias para la elaboración de la práctica. Entre las mismas, destacamos:
caret (Classification And REgresion Training): paquete para entrenar modelos de machine learning en R. Facilita, entre otros, el preprocesado de datos, la división de conjuntos de train/test, el ajuste de hiperparámetros…
tidyverse: colección de paquetes para ciencia de datos. Permite manipulación de datos, lectura de archivos y creación de gráficos.
ggplot2: creación de gráficos.
gridExtra: permite combinar varios ggplot en una sola figura.
También cargamos el dataset
y vemos su contenido. Observamos una tabla de 690 filas y 16 columnas.
Las columnas V1, V4, V5,
V6, V7, V9, V10,
V12, V13 y V16 contienen
caracteres, lo que nos indica que son categóricas. El resto de columnas
son numéricas.
if (!requireNamespace("caret", quietly = TRUE)) {
install.packages("caret")
}
if (!requireNamespace("randomForest", quietly = TRUE)) {
install.packages("randomForest")
}
if (!requireNamespace("doParallel", quietly = TRUE)) {
install.packages("doParallel")
}
if (!requireNamespace("tidyverse", quietly = TRUE)) {
install.packages("tidyverse")
}
if (!requireNamespace("ggplot2", quietly = TRUE)) {
install.packages("ggplot2")
}
if (!requireNamespace("gridExtra", quietly = TRUE)) {
install.packages("gridExtra")
}
if (!requireNamespace("GGally", quietly = TRUE)) {
install.packages("GGally")
}
if (!requireNamespace("rpart", quietly = TRUE)) {
install.packages("rpart")
}
if (!requireNamespace("rpart.plot", quietly = TRUE)) {
install.packages("rpart.plot")
}
if (!requireNamespace("gbm", quietly = TRUE)) {
install.packages("gbm")
}
if (!requireNamespace("nnet", quietly = TRUE)) {
install.packages("nnet")
}
if (!requireNamespace("tidyr", quietly = TRUE)) {
install.packages("tidyr")
}
if (!requireNamespace("dplyr", quietly = TRUE)) {
install.packages("dplyr")
}
library(nnet)
library(gbm)
library(caret)
library(tidyverse)
library(ggplot2)
library(gridExtra)
library(GGally)
library(tidyr)
library(dplyr)
#cargamos dataset
url <- "https://archive.ics.uci.edu/static/public/27/credit+approval.zip"
temp <- tempfile()
download.file(url, temp)
unzip(temp, exdir = "./credit")
unlink(temp)
credit <- read.table("./credit/crx.data", sep = ",", na.strings ="?")
summary(credit)
## V1 V2 V3 V4
## Length:690 Min. :13.75 Min. : 0.000 Length:690
## Class :character 1st Qu.:22.60 1st Qu.: 1.000 Class :character
## Mode :character Median :28.46 Median : 2.750 Mode :character
## Mean :31.57 Mean : 4.759
## 3rd Qu.:38.23 3rd Qu.: 7.207
## Max. :80.25 Max. :28.000
## NA's :12
## V5 V6 V7 V8
## Length:690 Length:690 Length:690 Min. : 0.000
## Class :character Class :character Class :character 1st Qu.: 0.165
## Mode :character Mode :character Mode :character Median : 1.000
## Mean : 2.223
## 3rd Qu.: 2.625
## Max. :28.500
##
## V9 V10 V11 V12
## Length:690 Length:690 Min. : 0.0 Length:690
## Class :character Class :character 1st Qu.: 0.0 Class :character
## Mode :character Mode :character Median : 0.0 Mode :character
## Mean : 2.4
## 3rd Qu.: 3.0
## Max. :67.0
##
## V13 V14 V15 V16
## Length:690 Min. : 0 Min. : 0.0 Length:690
## Class :character 1st Qu.: 75 1st Qu.: 0.0 Class :character
## Mode :character Median : 160 Median : 5.0 Mode :character
## Mean : 184 Mean : 1017.4
## 3rd Qu.: 276 3rd Qu.: 395.5
## Max. :2000 Max. :100000.0
## NA's :13
La primera labor que debemos llevar a cabo al empezar a trabajar con un dataset desconocido es entender la información que este contiene. En esta fase, conocida como Análisis Exploratorio de Datos, llevaremos a cabo las siguientes tareas:
Etiquetado de los atributos: en vista de la ausencia de etiquetas interpretables para los atributos del dataset, en esta primera etapa hemos buscado entender el significado de cada uno de ellos y renombrarlos de manera más intuitiva.
Análisis monovariable: a continuación, hemos analizado las distribuciones seguidas por cada uno de los atributos, estudiando también los valores máximos, mínimos, medias y cuartiles de aquellos campos numéricos.
Análisis multivariable: finalmente, hemos realizado un estudio de cómo distintos atributos se relacionan entre sí.
En un primer momento, intentamos averiguar los significados de los
campos de forma manual. Conseguimos deducir que el atributo
V2 representa la edad del solicitante, percatándonos de que
los decimales correspondían con múltiplos de 1/12 (por ejemplo, el valor
34.083 se corresponde con la edad de 34 años y 1 mes).
Otro campo que pudimos deducir fue el de Ingresos (V15),
basándonos en la distribución que tomaban sus valores (y que veremos más
adelante).
Finalmente, decidimos investigar en foros online para buscar el significado del resto de atributos. Ahí es cuando encontramos el siguiente proyecto de Kaggle.
Con esta información, renombramos las columnas:
colnames(credit)[colnames(credit) == "V1"] ="Sexo"
colnames(credit)[colnames(credit) == "V2"] ="Edad"
colnames(credit)[colnames(credit) == "V3"] ="Deuda"
colnames(credit)[colnames(credit) == "V4"] ="Estado_civil"
colnames(credit)[colnames(credit) == "V5"] ="Es_cliente"
colnames(credit)[colnames(credit) == "V6"] ="Nivel_educativo"
colnames(credit)[colnames(credit) == "V7"] ="Etnia"
colnames(credit)[colnames(credit) == "V8"] ="Anos_cotizados"
colnames(credit)[colnames(credit) == "V9"] ="Impago_previo"
colnames(credit)[colnames(credit) == "V10"] ="Trabaja"
colnames(credit)[colnames(credit) == "V11"] ="Calificacion_crediticia"
colnames(credit)[colnames(credit) == "V12"] ="Licencia_de_conducir"
colnames(credit)[colnames(credit) == "V13"] ="Ciudadano"
colnames(credit)[colnames(credit) == "V14"] ="Codigo_postal"
colnames(credit)[colnames(credit) == "V15"] ="Ingresos"
colnames(credit)[colnames(credit) == "V16"] ="Aprobado"
También distinguimos entre campos categóricos y numéricos. Los datos categóricos se deben tratar como el tipo factor. Algunos campos numéricos conviene visualizarlos en escala logarítmica.
campos = 1:15
campos_numericos = c(2, 3, 8, 11, 14, 15)
campos_log = c(8, 11, 15)
campos_no_log = setdiff(campos_numericos, campos_log)
campos_categoricos = setdiff(campos, campos_numericos)
credit[campos_categoricos] <- lapply(credit[campos_categoricos], FUN = factor)
lapply(credit[campos_categoricos], FUN = levels)
## $Sexo
## [1] "a" "b"
##
## $Estado_civil
## [1] "l" "u" "y"
##
## $Es_cliente
## [1] "g" "gg" "p"
##
## $Nivel_educativo
## [1] "aa" "c" "cc" "d" "e" "ff" "i" "j" "k" "m" "q" "r" "w" "x"
##
## $Etnia
## [1] "bb" "dd" "ff" "h" "j" "n" "o" "v" "z"
##
## $Impago_previo
## [1] "f" "t"
##
## $Trabaja
## [1] "f" "t"
##
## $Licencia_de_conducir
## [1] "f" "t"
##
## $Ciudadano
## [1] "g" "p" "s"
#pasamos el campo aprobado también a factor y renombramos sus valores
credit <- credit %>%
mutate(Aprobado = factor(Aprobado, levels = c("+", "-"), labels = c("Aprobado", "Rechazado")))
Añadimos los valores que faltan que pueden tomar algunas de las variables, en este caso a la variable estado civil le falta el valor “t”.
levels(credit$Estado_civil) <- c(levels(credit$Estado_civil), "t")
Para llevar a cabo el análisis monovariable, distinguiremos entre los atributos numéricos y categóricos, pues la forma de tratarlos es completamente distinta.
Para los atributos numéricos, es interesante como punto de partida mostrar su valor mínimo, máximo, los cuartiles, la media y la mediana.
#tabla con los atributos
atributosNum <- credit[,c(2,3,8,11,14,15)]
summary(atributosNum)
## Edad Deuda Anos_cotizados Calificacion_crediticia
## Min. :13.75 Min. : 0.000 Min. : 0.000 Min. : 0.0
## 1st Qu.:22.60 1st Qu.: 1.000 1st Qu.: 0.165 1st Qu.: 0.0
## Median :28.46 Median : 2.750 Median : 1.000 Median : 0.0
## Mean :31.57 Mean : 4.759 Mean : 2.223 Mean : 2.4
## 3rd Qu.:38.23 3rd Qu.: 7.207 3rd Qu.: 2.625 3rd Qu.: 3.0
## Max. :80.25 Max. :28.000 Max. :28.500 Max. :67.0
## NA's :12
## Codigo_postal Ingresos
## Min. : 0 Min. : 0.0
## 1st Qu.: 75 1st Qu.: 0.0
## Median : 160 Median : 5.0
## Mean : 184 Mean : 1017.4
## 3rd Qu.: 276 3rd Qu.: 395.5
## Max. :2000 Max. :100000.0
## NA's :13
Vamos a mostrar ahora gráficamente estos atributos, junto con un histograma que nos sirve como primera toma de contacto con la distribución que siguen dichos atributos.
plots <- list()
for (i in campos_numericos) {
var <- names(credit)[i]
var_data <- na.omit(credit[[var]])
if (length(var_data) < 2) next
df_temp <- data.frame(valor = var_data)
p <- ggplot(df_temp, aes(x = valor)) +
geom_histogram(bins = 30, fill = "skyblue", color = "black") +
ggtitle(var) +
theme_minimal()
plots[[var]] <- p
}
grid.arrange(grobs = plots, ncol = 3)
Claramente hay variables como la de Ingresos o
Calificacion_crediticia que pierden mucha información al
representarse con el histograma. Utilizando escala logarítmica se puede
distinguir mejor la población:
plots <- list()
for (i in campos_log) {
var <- names(credit)[i]
var_data <- na.omit(credit[[var]])
var_data_log <- var_data[var_data > 0]
if (length(var_data_log) < 2) next
df_temp <- data.frame(valor = log(var_data_log))
p <- ggplot(df_temp, aes(x = valor)) +
geom_histogram(bins = 30, fill = "skyblue", color = "black") +
ggtitle(paste0(var, " (log)")) +
theme_minimal()
plots[[var]] <- p
}
grid.arrange(grobs = plots, ncol = 3)
En el caso del atributo Calificacion_crediticia, vemos
que hay una gran cantidad de muestras que tienen el valor 0 en este
campo. Esto se puede deber a que el banco no conoce este dato o bien que
la mayor parte de la gente acaba de abrirse la cuenta del banco y de
momento no tiene una calificación crediticia positiva.
También vemos una cantidad elevada de personas con cero ingresos. Igual que en el caso anterior, seguramente se deba a desconocimiento del valor, pues no es frecuente que una cantidad tan elevada de la población tenga 0 ingresos.
Ahora queremos comprobar si alguno de los atributos sigue una distribución normal. Para ello utilizamos qqplot para comparar la función de distribución de nuestros atributos con la de la distribución normal.
qqplots <- list()
for (i in campos_numericos) {
var <- names(credit)[i]
var_data <- na.omit(credit[[i]])
if (length(var_data) < 3) next
df_temp <- data.frame(valor = var_data)
#Q-Q plot usando stat_qq
p <- ggplot(df_temp, aes(sample = valor)) +
stat_qq() +
stat_qq_line(color = "red") +
ggtitle(paste(var)) +
theme_minimal()
qqplots[[var]] <- p
}
#mostramos todos los Q-Q plots en una cuadrícula
grid.arrange(grobs = qqplots, ncol = 3)
Llegamos a la conclusión de que ninguno de los atributos sigue una distribución normal. De seguir esta distribución, deberíamos ver una coincidencia entre los cuartiles reales (puntos) y la línea de cuartiles (roja) para una distribución normal de misma media y desviación típica.
La que más se acerca sería la del código postal. Esto podría tener cierto sentido debido a que la mayor parte de los solicitantes vivirán cerca del banco en cuestión. Teniendo en cuenta que lugares cercanos suelen tener códigos postales próximos, tiene mucho sentido que cuanto más nos alejamos del banco (y por consecuente de su código postal), menos solicitantes habrá.
Siguiendo el razonamiento anterior, parece convenienente comprobar si las variables que anteriormente visualizamos en escala logarítmica, siguen una distribución lognormal. Para ello simplemente repetimos el procedimiento anterior con las variables tras aplicarle el logaritmo:
qqplots <- list()
for (i in campos_log) {
var <- names(credit)[i]
var_data <- na.omit(credit[[var]])
var_data_log <- var_data[var_data > 0]
if (length(var_data_log) < 2) next
df_temp <- data.frame(valor = log(var_data_log))
#Q-Q plot usando stat_qq
p <- ggplot(df_temp, aes(sample = valor)) +
stat_qq() +
stat_qq_line(color = "red") +
ggtitle(paste(var)) +
theme_minimal()
qqplots[[var]] <- p
}
#mostramos todos los Q-Q plots en una cuadrícula
grid.arrange(grobs = qqplots, ncol = 3)
Vemos que estas tres variables se ajustan mucho mejor a una
distribución lognormal que a una normal. Ahora vamos a centrarnos en dos
predictores en concreto: Ingresosy Deuda.
Ingresos#eliminamos los valores NA
ingresos <- na.omit(credit$Ingresos)
#Visualizamos los estadísticos básicos
summary(ingresos)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 0.0 5.0 1017.4 395.5 100000.0
sd(ingresos)
## [1] 5210.103
# Histograma y Boxplot
p1 <- ggplot(data.frame(Ingresos = ingresos), aes(x = Ingresos)) +
geom_histogram(bins = 30, fill = "steelblue", color = "black") +
labs(title = "Histograma de Ingresos")
p2 <- ggplot(data.frame(Ingresos = ingresos), aes(y = Ingresos)) +
geom_boxplot(fill = "orange") +
labs(title = "Boxplot de Ingresos")
gridExtra::grid.arrange(p1, p2, ncol = 2)
La variable Ingresos presenta una distribución altamente
asimétrica, con una gran concentración de valores bajos y una cola larga
hacia la derecha. Esto se confirma con el histograma y el boxplot, donde
se observan muchos valores en cero y varios outliers.
Para comprobar si sigue una distribución normal podemos utilizar el test de Shapiro-Wilk:
# Test de Shapiro-Wilk para normalidad
shapiro.test(ingresos)
##
## Shapiro-Wilk normality test
##
## data: ingresos
## W = 0.16985, p-value < 2.2e-16
# Log-transformación (se suman +1 para evitar log(0))
log_ingresos <- log(ingresos + 1)
# Test de Shapiro-Wilk sobre log(Ingresos)
shapiro.test(log_ingresos)
##
## Shapiro-Wilk normality test
##
## data: log_ingresos
## W = 0.82425, p-value < 2.2e-16
La prueba de Shapiro-Wilk sobre los ingresos originales da como resultado un p-valor muy bajo, indicando que no siguen una distribución normal.
Tras aplicar una transformación logarítmica (log(Ingresos + 1)), la distribución mejora notablemente en forma, aunque el p-valor obtenido sigue siendo muy bajo, lo que nos confirma que tampoco sigue una distribución lognormal. El Q-Q plot también muestra una mejor alineación con la recta teórica tras la transformación.
# Q-Q plot para Ingresos y log(Ingresos)
par(mfrow = c(1,2))
qqnorm(ingresos, main = "Q-Q plot Ingresos")
qqline(ingresos, col = "red")
qqnorm(log_ingresos, main = "Q-Q plot log(Ingresos)")
qqline(log_ingresos, col = "red")
Nuestra siguiente intención es comprobar si este predictor sirve como una buena primera aproximación para distinguir entre personas a las que se le concedió el crédito y a las que no. Para ello, dividimos la población total entre aprobados y no aprobados y representamos con un boxplot cada una.
# Eliminar NA
df_ingresos <- credit[!is.na(credit$Ingresos) & !is.na(credit$Aprobado), ]
# Boxplot de Ingresos por Clase
ggplot(df_ingresos, aes(x = Aprobado, y = Ingresos, fill = Aprobado)) +
geom_boxplot() +
labs(title = "Distribución de Ingresos según Aprobado",
x = "Aprobado (Crédito aprobado o no)",
y = "Ingresos") +
theme_minimal()
Al tratarse la mayoría de valores muy cercanos a cero perdemos casi toda información al visualizar la gráfica. No obstante sí que nos damos cuenta de que a todos los outliers (personas con ingresos muy por encima de la media) se les concede el crédito, lo que nos hace pensar que puede ser un buen predictor. Para visualizar mejor el boxplot, transformamos los datos a escala logaritmica:
df_ingresos$log_ingresos <- log(df_ingresos$Ingresos + 1)
ggplot(df_ingresos, aes(x = Aprobado, y = log_ingresos, fill = Aprobado)) +
geom_boxplot() +
labs(title = "Distribución log(Ingresos + 1) según Aprobado",
x = "Aprobado",
y = "log(Ingresos + 1)") +
theme_minimal()
En este diagrama sí que se aprecia mucho mejor que a la gente que se le concedió el crédito tiene de media muchos más ingresos, confirmando así nuestra teoría. Podemos mostrar también las funciones de densidad para cada una de las clases, en escala logarítmica para una mejor visualización:
ggplot(df_ingresos, aes(x = log_ingresos, fill = Aprobado)) +
geom_density(alpha = 0.5) +
labs(title = "Distribución de log(Ingresos+1) según Aprobado",
x = "log(Ingresos+1)",
y = "Densidad") +
theme_minimal()
Deuda# Eliminamos NA
deuda <- na.omit(credit$Deuda)
# Estadísticos descriptivos
summary(deuda)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.000 2.750 4.759 7.207 28.000
sd(deuda)
## [1] 4.978163
# Histograma y Boxplot
p1 <- ggplot(data.frame(Deuda = deuda), aes(x = Deuda)) +
geom_histogram(bins = 30, fill = "skyblue", color = "black") +
labs(title = "Histograma de Deuda")
p2 <- ggplot(data.frame(Deuda = deuda), aes(y = Deuda)) +
geom_boxplot(fill = "skyblue") +
labs(title = "Boxplot de Deuda")
gridExtra::grid.arrange(p1, p2, ncol = 2)
La variable Deuda presenta una distribución asimétrica,
aunque menos extrema que Ingresos. El histograma muestra
una acumulación de valores bajos y una cola hacia la derecha. El boxplot
sugiere la presencia de varios valores atípicos por encima del tercer
cuartil.
# Prueba de normalidad
shapiro.test(deuda)
##
## Shapiro-Wilk normality test
##
## data: deuda
## W = 0.83025, p-value < 2.2e-16
# Log-transformación si hay valores positivos
log_deuda <- log(deuda + 1)
shapiro.test(log_deuda)
##
## Shapiro-Wilk normality test
##
## data: log_deuda
## W = 0.96347, p-value = 4.345e-12
Ambas pruebas de Shapiro-Wilk indican que ni Deuda ni
log(Deuda + 1) siguen una distribución normal. No obstante, el logaritmo
ayuda a mejorar la simetría:
# Q-Q plots
par(mfrow = c(1,2))
qqnorm(deuda, main = "Q-Q plot Deuda")
qqline(deuda, col = "red")
qqnorm(log_deuda, main = "Q-Q plot log(Deuda)")
qqline(log_deuda, col = "red")
Queremos ahora ver si Deuda permite distinguir entre
personas a las que se les concedió el crédito y a las que no:
df_deuda <- credit[!is.na(credit$Deuda) & !is.na(credit$Aprobado), ]
df_deuda$log_deuda <- log(df_deuda$Deuda + 1)
ggplot(df_deuda, aes(x = Aprobado, y = log_deuda, fill = Aprobado)) +
geom_boxplot() +
labs(title = "Distribución log(Deuda + 1) según Aprobado",
x = "Aprobado",
y = "log(Deuda + 1)") +
theme_minimal()
El boxplot log-transformado permite ver que los casos aprobados tienden a tener mayor deuda, lo cual es contraintuitivo, pero puede estar relacionado con el hecho de que esas personas también tienen mayores ingresos y capacidad de pago. Este comportamiento debe analizarse junto con otros predictores.
ggplot(df_deuda, aes(x = Deuda, fill = Aprobado)) +
geom_density(alpha = 0.5) +
labs(title = "Distribución de Deuda según Aprobado",
x = "Deuda",
y = "Densidad") +
theme_minimal()
Este gráfico muestra que, aunque las distribuciones se solapan, el grupo aprobado presenta una mayor densidad en valores altos de deuda, reforzando la idea de que tener deuda no impide ser aprobado, siempre que venga acompañada de otros factores como ingresos estables.
Podemos verlo en escala logarítmica de manera más evidente:
ggplot(df_deuda, aes(x = log_deuda, fill = Aprobado)) +
geom_density(alpha = 0.5) +
labs(title = "Distribución de Deuda según Aprobado",
x = "Deuda",
y = "Densidad") +
theme_minimal()
El estudio de campos categóricos es limitado en comparación con el de campos numéricos. Al no existir orden ni distancia entre valores de los atributos, perdemos la ayuda de la mayor parte de herramientas estadísticas. Sin embargo, esto no significa que no podamos extraer información útil de estos campos. Podemos comenzar nuetro análisis graficando un histograma de cada variable categórica.
summary(credit[campos_categoricos])
## Sexo Estado_civil Es_cliente Nivel_educativo Etnia Impago_previo
## a :210 l : 2 g :519 c :137 v :399 f:329
## b :468 u :519 gg : 2 q : 78 h :138 t:361
## NA's: 12 y :163 p :163 w : 64 bb : 59
## t : 0 NA's: 6 i : 59 ff : 57
## NA's: 6 aa : 54 j : 8
## (Other):289 (Other): 20
## NA's : 9 NA's : 9
## Trabaja Licencia_de_conducir Ciudadano
## f:395 f:374 g:625
## t:295 t:316 p: 8
## s: 57
##
##
##
##
categorical_vars <- names(credit)[campos_categoricos]
# Create a plotting function for categorical variables
plot_categorical <- function(var_name) {
ggplot(credit, aes(x = .data[[var_name]])) +
geom_bar(fill = "skyblue", color = "black") +
labs(title = paste("Distribución de", var_name),
x = var_name, y = "Cantidad") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
}
categorical_plots <- map(categorical_vars, plot_categorical)
grid.arrange(grobs = categorical_plots, ncol = 2)
Nos gustaría también ver cómo los valores que toman estos campos
afectan a la variable a predecir. Para ello podemos hacer uso de la
función spineplot, que nos muestra para cada valor de una
variable tanto su poblacion (ancho) como su tasa de aprobados (alto de
la barra azul).
par(mfrow = c(3,3)) #plots en cuadrícula 3x3
for (var in categorical_vars) {
spineplot(x = credit[[var]],
y = credit$Aprobado,
ylevels = c("Rechazado", "Aprobado"),
xlab = var,
ylab = "Resultado",
main = paste("Aprobados por", var),
col = c("#46cac1", "#ff8080"),
border = c("#26958d", "#b13636"))
}
par(mfrow = c(1, 1))
Como punto de interés adicional, representemos también etnia frente a ciudadanía:
spineplot(credit$Etnia, credit$Ciudadano,
xlab = "Etnia",
ylab = "Es ciudadano",
main = "Ciudadanía por Etnia",
col = c("#aeaeae", "#848484", "#575757"))
Con los datos representados podemos realizar varias deducciones sobre el contexto del dataset y los significados de las etiquetas de algunos campos. Por ejemplo:
En el plot de Etnia frente a Ciudadanía,
vemos que no existen grandes variaciones en las distribuciones de
ciudadanía entre distintas etnias. Esto nos hace pensar que no estamos
tratando con datos provenientes de un único país. Si fuese así, la etnia
mayoritaria tendría probablemente un porcentaje de miembros ciudadanos
muy superior al resto. Debemos suponer entonces, que aunque los datos
provengan de un banco japonés, sus clientes están distribuidos por
varios países, y el campo Es_ciudadano hace referencia a la
cuidadanía respecto al país en el que residen.
Trabaja: puesto que una persona con trabajo tiene mucha
más facilidad para recibir aprobación por parte de un banco, resulta
bastante claro que “t” y “f” se corresponden con tener y no tener
trabajo, respectivamente.
Impago previo: vemos de manera bastante inmediata que el
valor “f” representa presencia de impagos por su terrible impacto en la
tasa de aprobado (al contrario de lo que el nombre del campo pueda dar a
entender). “t” significaría entonces ausencia de impagos. El valor de
esta variable es un factor muy importante a la hora de decidir si se
aprueba o no la solicitud. Tanto así que un predictor que utilice
simplemente esta variable para clasificar, obtendría una tasa de
aciertos de (306+284)/690 ~= 0.855, bastante cercano a lo que
obtendremos con técnicas más avanzadas.
tab <- table(credit$Impago_previo, credit$Aprobado)
df_plot <- as.data.frame(tab)
names(df_plot) <- c("Impago_previo", "Aprobado", "Count")
ggplot(df_plot, aes(x = Impago_previo, y = Aprobado, fill = Count)) +
geom_tile(color = "white") +
geom_text(aes(label=Count), color = "black", size = 6) +
scale_fill_gradient(low = "white", high = "steelblue") +
labs(title = "Tabla de frecuencias de impago previo frente a aprobado",
x = "Impago previo",
y = "Aprobado") +
theme_minimal(base_size = 12)
Analizar todas las posibles agrupaciones de atributos sería una tarea inabarcable. Al contar con 16 de los mismos, existen un total de 216 = 65536 combinaciones de conjuntos. Es por ello que incluimos solo análisis de grupos de atributos que consideramos que pueden tener importancia.
Los objetivos del análisis multivariable es detectar patrones complejos no visibles en análisis univariantes. Nos puede servir posteriormente para reducir dimensiones o justificar otras decisiones de preprocesado.
Utilizando un pairwise correlation plot con la
función del paquete GGally vamos a visualizar las relaciones entre todos
los pares de variables numéricas. Distinguimos entre créditos aprobados
(rojo) y denegados (azul). No incluimos el Código_postal ya
que no es un valor que mantenga un concepto de escala.
Dado que las variables Deuda,
Anos_cotizados y Calificacion_crediticia
presentan distribuciones muy sesgadas hacia valores bajos, aplicamos una
transformación logarítmica. Así, se estabilizamos las varianzas y
favorecemos el análisis de regresiones lineales.
#función para mostrar R^2 en los paneles superiores
cor_plus_r2_fun <- function(data, mapping, ...) {
x <- eval_data_col(data, mapping$x)
y <- eval_data_col(data, mapping$y)
class_var <- data$Aprobado # Asumimos que la clase es esta
df <- data.frame(x = x, y = y, clase = class_var)
df <- df[complete.cases(df), ]
#Correlaciones
r_total <- cor(df$x, df$y)
r_plus <- tryCatch(cor(df$x[df$clase == "Aprobado"], df$y[df$clase == "Aprobado"]), error = function(e) NA)
r_minus <- tryCatch(cor(df$x[df$clase == "Rechazado"], df$y[df$clase == "Rechazado"]), error = function(e) NA)
#R2
r2 <- tryCatch(summary(lm(y ~ x, data = df))$r.squared, error = function(e) NA)
#Formato
txt <- paste0(
"Corr: ", sprintf("%.3f", r_total), "\n",
"-: ", sprintf("%.3f", r_minus), "\n",
"+: ", sprintf("%.3f", r_plus), "\n",
"R²: ", sprintf("%.3f", r2)
)
ggplot(data = df, aes(x = x, y = y)) +
annotate("text", x = median(df$x, na.rm = TRUE), y = median(df$y, na.rm = TRUE),
label = txt, size = 3.5, hjust = 0.5) +
theme_void()
}
#Subconjunto con las variables numéricas
subset_credit <- credit[, c("Edad", "Deuda", "Anos_cotizados", "Calificacion_crediticia", "Ingresos", "Aprobado")]
subset_credit$Deuda <- log10(subset_credit$Deuda + 1)
subset_credit$Anos_cotizados <- log10(subset_credit$Anos_cotizados + 1)
subset_credit$Ingresos <- log10(subset_credit$Ingresos + 1)
subset_credit$Calificacion_crediticia <- log10(subset_credit$Calificacion_crediticia + 1)
#Graficamos con color por clase
ggpairs(subset_credit,
mapping = aes(color = Aprobado, alpha = 0.5),
columns = 1:5,
upper = list(continuous = cor_plus_r2_fun),
lower = list(continuous = wrap("points", alpha = 0.4)),
title = "Relaciones entre atributos numéricos")
## Warning: Removed 12 rows containing non-finite outside the scale range
## (`stat_density()`).
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).
Como puede observarse en la matriz de correlaciones, la diagonal corresponde a las distribuciones marginales que mostrábamos en el análisis monovariable. El triángulo inferior contiene los diagramas de dispersión para cada par de variables y clase (aprobado o denegado). Por último, el triángulo superior derecho contiene los coeficientes de correlación de Pearson, tanto de las distribuciones globales como las separadas por clase.
Las correlaciones son, en su mayoría, positivas y de pequeña
magnitud. Tiene sentido, por ejemplo, que a mayor edad tenga una
persona, más larga sea su vida laboral. Estas correlaciones son más
intensas en la clase +. Esto sugiere que los casos
aprobados siguen patrones más regulares. En el caso de la relación
Edad-Anos_cotizados, podríamos interpretar que
a la hora de la concesión de un crédito a una persona joven, se es mas
laxo con una vida laboral más corta. No obstante, si nos fijamos en los
gráficos de dispersión y en el R2, podemos concluir que los
datos no son los idóneos para una regresión lineal.
Como separador, el par Edad-Deuda es el
peor de todos, ya que las distribuciones de ambas clases se solapan
demasiado. Todos los pares en los que aparece la calificación crediticia
son buenas elecciones de separadores; no por ser un buen par, sino
porque la calificación crediticia es en sí misma un buen separador. El
par Deuda-Anos_cotizados es quizá el mejor
separador bidimensional de la tabla. Se puede apreciar que un corte
diagonal (recta y=1.5-x) se consigue una buena separación de los puntos
rojos y azules.
subset_credit$Esfuerzo_financiero <- subset_credit$Deuda + subset_credit$Anos_cotizados + abs(subset_credit$Ingresos-1)*0.375
ggplot(subset_credit, aes(x = Esfuerzo_financiero, fill = Aprobado)) +
geom_density(alpha = 0.5) +
theme_minimal() +
labs(title = "Distribución del índice combinado log10(Deuda+1) + log10(Años cotizados+1) + |Igresos-1|*0.375" ,
x = "Índice esfuerzo financiero", y = "Densidad")
Hemos llamado a la suma de ambos atributos (en escala logarítmica)
mas el término |Ingresos-1|*0.375 Esfuerzo_financiero. La
distribución marginal de este nuevo valor tiende a valores cercanos a
cero en el caso de los créditos denegados y a una distribución simétrica
centrada a la mitad del eje. A pesar del solapamiento, el patrón sugiere
que este nuevo atributo podría tener valor como variable predictiva.
Podría ser interesante como incorporación a los modelos de
aprendizaje.
La calificación crediticia es, a priori, el mejor clasificador numérico monovariable. No obstante, sigue habiendo unos pocos casos en los que una calificación alta no implica la concesión del crédito. En esta subsección vamos a intentar encontrar cuáles son las variables que determinan la concesión del crédito en los casos de un índice crediticio alto.
Para empezar, filtramos el dataset para contener únicamente información sobre los créditos en los que el solicitante tiene un índice crediticio superior a 3.
credit_alto <- subset(credit, Calificacion_crediticia > 3)
summary(credit_alto)
## Sexo Edad Deuda Estado_civil Es_cliente
## a:54 Min. :16.08 Min. : 0.000 l: 0 g :133
## b:97 1st Qu.:24.00 1st Qu.: 1.750 u:133 gg: 0
## Median :33.17 Median : 5.665 y: 18 p : 18
## Mean :35.43 Mean : 6.739 t: 0
## 3rd Qu.:43.16 3rd Qu.:10.312
## Max. :68.67 Max. :28.000
##
## Nivel_educativo Etnia Anos_cotizados Impago_previo Trabaja
## c :34 v :80 Min. : 0.000 f: 10 f: 0
## q :26 h :40 1st Qu.: 1.000 t:141 t:151
## cc :14 bb :15 Median : 2.375
## w :13 ff : 9 Mean : 4.010
## x :13 z : 4 3rd Qu.: 5.250
## i : 9 dd : 1 Max. :28.500
## (Other):42 (Other): 2
## Calificacion_crediticia Licencia_de_conducir Ciudadano Codigo_postal
## Min. : 4.000 f:82 g:149 Min. : 0.0
## 1st Qu.: 6.000 t:69 p: 0 1st Qu.: 0.0
## Median : 8.000 s: 2 Median :100.0
## Mean : 9.344 Mean :141.4
## 3rd Qu.:11.000 3rd Qu.:220.0
## Max. :67.000 Max. :583.0
## NA's :2
## Ingresos Aprobado
## Min. : 0.0 Aprobado :135
## 1st Qu.: 7.5 Rechazado: 16
## Median : 540.0
## Mean : 1919.4
## 3rd Qu.: 1802.5
## Max. :51100.0
##
El resultado es una tabla de 151 entradas. El perfil medio de estas solicitudes, más allá de la calificación crediticia, se deferencia del perfil general en cuatro ámbitos:
Todos tienen un puestro de trabajo.
La deuda es mayor (6.739 vs 4.759 de media).
Más años cotizados (1 en el primer cuartil vs 0.165).
Más ingresos (540 vs 5 de media).
Si repetimos la matriz de correlaciones para este dataset:
#función para mostrar R^2 en los paneles superiores
cor_plus_r2_fun <- function(data, mapping, ...) {
x <- eval_data_col(data, mapping$x)
y <- eval_data_col(data, mapping$y)
class_var <- data$Aprobado # Asumimos que la clase es esta
df <- data.frame(x = x, y = y, clase = class_var)
df <- df[complete.cases(df), ]
#Correlaciones
r_total <- cor(df$x, df$y)
r_plus <- tryCatch(cor(df$x[df$clase == "Aprobado"], df$y[df$clase == "Aprobado"]), error = function(e) NA)
r_minus <- tryCatch(cor(df$x[df$clase == "Rechazado"], df$y[df$clase == "Rechazado"]), error = function(e) NA)
#R^2 global
r2 <- tryCatch(summary(lm(y ~ x, data = df))$r.squared, error = function(e) NA)
#Formato
txt <- paste0(
"Corr: ", sprintf("%.3f", r_total), "\n",
"-: ", sprintf("%.3f", r_minus), "\n",
"+: ", sprintf("%.3f", r_plus), "\n",
"R²: ", sprintf("%.3f", r2)
)
ggplot(data = df, aes(x = x, y = y)) +
annotate("text", x = median(df$x, na.rm = TRUE), y = median(df$y, na.rm = TRUE),
label = txt, size = 3.5, hjust = 0.5) +
theme_void()
}
# Subconjunto con las variables numéricas
subset_credit <- credit_alto[, c("Edad", "Deuda", "Anos_cotizados", "Calificacion_crediticia", "Ingresos", "Aprobado")]
subset_credit$Deuda <- log10(subset_credit$Deuda + 1)
subset_credit$Anos_cotizados <- log10(subset_credit$Anos_cotizados + 1)
subset_credit$Ingresos <- log10(subset_credit$Ingresos + 1)
subset_credit$Calificacion_crediticia <- log10(subset_credit$Calificacion_crediticia + 1)
#graficamos con color por clase
ggpairs(subset_credit,
mapping = aes(color = Aprobado, alpha = 0.5),
columns = 1:5,
upper = list(continuous = cor_plus_r2_fun),
lower = list(continuous = wrap("points", alpha = 0.4)),
title = "Relaciones entre atributos numéricos con alta Calificacion_crediticia")
los resultados obtenidos sugieren que los ingresos del solicitante
son un factor clave en la concesión o no del crédito, resultando en una
denegación aquellos perfiles con unos ingresos cercanos a 10 unidades.
También destaca la diferencia entre una deuda alta y baja. De hecho, si
utilizamos el Esfuerzo_financiero que hemos definido antes
obtenemos el siguiente resultado:
subset_credit$Esfuerzo_financiero <- abs(subset_credit$Ingresos-1)*0.375 + subset_credit$Anos_cotizados + subset_credit$Deuda
ggplot(subset_credit, aes(x = Esfuerzo_financiero, fill = Aprobado)) +
geom_density(alpha = 0.5) +
theme_minimal() +
labs(title = "Distribución log10(Deuda+1) + log10(Años cotizados+1) + |Ingresos-1|*0.375",
x = "Índice esfuerzo financiero", y = "Densidad")
La separación es aún mejor que en el caso general. Podemos ver el rendimiento de este predictor para este subdataset:
subset_credit$Prediccion_simple <- ifelse(subset_credit$Esfuerzo_financiero > 1.6, "Aprobado", "Rechazado")
table(Real = subset_credit$Aprobado, Predicho = subset_credit$Prediccion_simple)
## Predicho
## Real Aprobado Rechazado
## Aprobado 102 33
## Rechazado 3 13
#creamos tabla de confusión en formato largo
conf_mat <- table(Real = subset_credit$Aprobado, Predicho = subset_credit$Prediccion_simple)
conf_df <- as.data.frame(conf_mat)
#añadimos etiquetas
conf_df <- conf_df %>%
mutate(Tipo = case_when(
Real == "Aprobado" & Predicho == "Aprobado" ~ "Verdadero Positivo",
Real == "Rechazado" & Predicho == "Aprobado" ~ "Falso Positivo",
Real == "Rechazado" & Predicho == "Rechazado" ~ "Verdadero Negativo",
Real == "Aprobado" & Predicho == "Rechazado" ~ "Falso Negativo"
))
# Gráfico
ggplot(conf_df, aes(x = Predicho, y = Real, fill = Freq)) +
geom_tile(color = "white") +
geom_text(aes(label = paste0(Freq, "\n", Tipo)), size = 4) +
scale_fill_gradient(low = "white", high = "steelblue") +
labs(title = "Matriz de confusión para predictor basado en esfuerzo financiero",
x = "Predicción", y = "Clase real", fill = "Frecuencia") +
theme_minimal()
La matriz de confusión nos indica que un predictor que utilizase este nuevo campo como clasificador en los casos con un índice crediticio alto, se obtendría una tasa de acierto del 75%.
credit_alto$Prediccion_simple <- ifelse(as.character(credit_alto$Impago_previo) == "t", "Aprobado", "Rechazado")
table(Real = credit_alto$Aprobado, Predicho = credit_alto$Prediccion_simple)
## Predicho
## Real Aprobado Rechazado
## Aprobado 135 0
## Rechazado 6 10
# Crear tabla de confusión en formato largo
conf_mat <- table(Real = credit_alto$Aprobado, Predicho = credit_alto$Prediccion_simple)
conf_df <- as.data.frame(conf_mat)
# Añadir etiquetas
conf_df <- conf_df %>%
mutate(Tipo = case_when(
Real == "Aprobado" & Predicho == "Aprobado" ~ "Verdadero Positivo",
Real == "Rechazado" & Predicho == "Aprobado" ~ "Falso Positivo",
Real == "Rechazado" & Predicho == "Rechazado" ~ "Verdadero Negativo",
Real == "Aprobado" & Predicho == "Rechazado" ~ "Falso Negativo"
))
# Gráfico
ggplot(conf_df, aes(x = Predicho, y = Real, fill = Freq)) +
geom_tile(color = "white") +
geom_text(aes(label = paste0(Freq, "\n", Tipo)), size = 4) +
scale_fill_gradient(low = "white", high = "steelblue") +
labs(title = "Matriz de confusión para predictor basado en impago previo",
x = "Predicción", y = "Clase real", fill = "Frecuencia") +
theme_minimal()
En cambio, si utilizamos aquí el predictor
Impago_previo, vemos como el resultado aumenta al 96%.
Realizamos ya la división entre Train y Test:
#cargamos índices de train
credit.trainIdx <- readRDS("credit.trainIdx")
#separamos los dos conjuntos
credit.Datos.Train <- credit[credit.trainIdx, ]
credit.Datos.Test <- credit[-credit.trainIdx, ]
En primer lugar, veamos si hay valores nulos y cuántos hay en cada columna:
#Contamos los nulos con una sola línea
colSums(is.na(credit))
## Sexo Edad Deuda
## 12 12 0
## Estado_civil Es_cliente Nivel_educativo
## 6 6 9
## Etnia Anos_cotizados Impago_previo
## 9 0 0
## Trabaja Calificacion_crediticia Licencia_de_conducir
## 0 0 0
## Ciudadano Codigo_postal Ingresos
## 0 13 0
## Aprobado
## 0
En una gráfica podemos verlos mejor:
#vemos losnulos con gráfico de barras
library(ggplot2)
nulos <- colSums(is.na(credit))
nulos_df <- data.frame(Variable = names(nulos), Nulos = nulos)
nulos_df <- nulos_df[nulos_df$Nulos > 0, ] # Solo mostramos las columnas que tengan nulos
ggplot(nulos_df, aes(x = reorder(Variable, -Nulos), y = Nulos)) +
geom_bar(stat = "identity", fill = "tomato") +
labs(title = "Número de valores nulos por variable",
x = "Variable", y = "Cantidad de NA") +
theme_minimal() +
coord_flip()
Vemos que solo 7 de los atributos tienen valores nulos y además la cantidad de estos es muy pequeña (cercano al 1% en todos los casos). Es interesante saber cómo se distribuyen los valores nulos a lo largo de las muestras. Para ello vamos a ver la frecuencia del número de atributos faltantes. Ya trabajamos con los datos de Train, pues no debemos modificar los datos de Test a partir de ahora.
#numero de valores NA por fila
na_por_fila <- rowSums(is.na(credit.Datos.Train ))
#observaciones con al menos un NA
con_na <- na_por_fila[na_por_fila > 0]
length(con_na)
## [1] 37
table(con_na)
## con_na
## 1 2 3 5
## 27 2 2 6
library(ggplot2)
ggplot(data.frame(Faltantes = con_na), aes(x = Faltantes)) +
geom_bar(fill = "steelblue") +
labs(title = "Nº de atributos faltantes por observación",
x = "Cantidad de NA",
y = "Frecuencia") +
theme_minimal()
Vemos que a la mayoría de muestras con valores NA solo les falta un valor. Sin embargo, también es importante destacar que hay 6 muestras a las cuales les faltan 5 valores.
Como estrategia de limpieza, se han eliminado las observaciones del conjunto de entrenamiento que contenían más de 4 valores nulos. Esta decisión se justifica por el hecho de que esas observaciones presentan información muy incompleta, y su imputación generaría más ruido que beneficio. El resto de las observaciones, con pocos valores faltantes, se conservarán y se decidirá si se inputan o no posteriormente.
#eliminamos las observaciones con más de 4 atributos faltantes
credit.Datos.Train <- credit.Datos.Train[na_por_fila <= 4, ]
#verificamos tamaño tras limpieza
nrow(credit.Datos.Train)
## [1] 547
Los valores nulos deben tratarse adecuadamente antes de entrenar modelos de aprendizaje automático. Las posibles estrategias incluyen:
Eliminación de filas: útil si el número de nulos es pequeño y su eliminación no reduce significativamente el tamaño del dataset.
Imputación por la media/mediana/moda: recomendable para variables numéricas si el porcentaje de nulos es moderado.
Imputación por valor más frecuente: útil en variables categóricas con pocos niveles.
Modelado de imputación: usar modelos predictivos (por ejemplo mice, missForest) para estimar los valores faltantes, aunque es más costoso computacionalmente.
En nuestro caso, hemos decidido simplemente eliminar las filas con datos faltantes, pues como hemos visto en la gráfica anterior, hay muy pocos valores nulos en comparación con el número total de muestras que disponemos.
Utilizando facetas y los diagramas boxplot podemos saber de la existencia de outliers para todas las variables numéricas:
#Selecciamos la variables numéricas
datos_numericos <- credit.Datos.Train %>%
select(where(is.numeric)) %>%
mutate(fila = row_number())
datos_long <- pivot_longer(datos_numericos,
cols = -fila,
names_to = "Variable",
values_to = "Valor")
#Boxplots con facet para cada variable
ggplot(datos_long, aes(x = "", y = Valor)) +
geom_boxplot(outlier.color = "red", fill = "lightblue") +
facet_wrap(~ Variable, scales = "free", ncol = 3) +
labs(title = "Detección de outliers por variable numérica",
y = "Valor",
x = "") +
theme_minimal()
Vemos que la mayoría de atributos sí que presentan outliers. Sin embargo, no hemos encontrado diferencias notables entre descartarlos y mantenerlos, pues en la transformación de datos aplicamos la transformación de Yeo-Johnson que prácticamente acaba con todos los outliers al normalizar estas variables.
EsfuerzoFinancieroComo ya hemos comentado, creemos que el atributo EsfuerzoFinanciero puede ser útil pues da buenos resultados como predictor. Es por eso que finalmente decidimos añadirlo como variable para entrenar los modelos:
credit.Datos.Train$Esfuerzo_financiero <- log10(credit.Datos.Train$Deuda + 1) + log10(credit.Datos.Train$Anos_cotizados + 1) + abs(log10(credit.Datos.Train$Ingresos + 1)-1)*0.375
credit.Datos.Test$Esfuerzo_financiero <- log10(credit.Datos.Test$Deuda + 1) + log10(credit.Datos.Test$Anos_cotizados + 1) + abs(log10(credit.Datos.Test$Ingresos + 1)-1)*0.375
if (!"Esfuerzo_financiero" %in% names(credit.Datos.Train)[campos_numericos]) {
campos_numericos <- c(campos_numericos, which(names(credit.Datos.Train) == "Esfuerzo_financiero"))
}
A la hora de transformar los datos hay varias opciones posibles, las cuales hemos comparado para ver los resultados que obteníamos con cada una de ellas:
Sin transformaciones: Es la opción más sencilla, que consiste en utilizar las variables directamente como vienen de serie.
Transformación logarítmica: Esta transformación es muy útil para variables con una distribución asimétrica positiva (sesgada a la derecha), como ingresos, deuda o puntuaciones de crédito, donde la mayoría de los valores son pequeños pero hay algunos muy grandes.
BoxCox: Intenta encontrar la mejor transformación para que una variable siga una distribución normal. La principal limitación es que no acepta valores negativos.
Yeo-Johnson: Es una extensión de la transformación Box-Cox que sí permite trabajar con valores negativos y con valores iguales a 0.
Tras probar las 4 alternativas, nos dimos cuenta que con Yeo-Johnson se obtenían ligeras mejoras sobre todo en el modelo gbm. El resto obtenían resultados muy similares con todas las transformaciones, por lo que decidimos quedarnos trabajar con la transformación de Yeo-Johnson para no tener que distinguir casos para cada uno de los modelos:
vars_to_transform <- names(credit.Datos.Train)[campos_numericos]
for (var in vars_to_transform) {
dx <- credit.Datos.Train[, var, drop = FALSE]
dx
dxTest <- credit.Datos.Test[, var, drop = FALSE]
#aplicamos la transformación Yeo-Johnson
preProcValues <- preProcess(dx, method = "YeoJohnson")
transformed <- predict(preProcValues, dx)
transformedTest <- predict(preProcValues, dxTest)
#sustituimos la columna original por la transformada
credit.Datos.Train[, var] <- transformed[[var]]
credit.Datos.Test[, var] <- transformedTest[[var]]
}
rm(dx, transformed, transformedTest, dxTest, var, vars_to_transform, preProcValues)
PCA es una técnica de análisis de datos utilizada para representar un conjunto de datos multidimensionales mediante un conjunto más reducido de variables generadas como combinaciones lineales de sus componentes. Esto se consigue buscando la proyección según la cual los datos presentan una mayor varianza. En nuestro caso solo 5 de los campos son numéricos. Los campos que antes representábamos en escala logarítmica se pasan también en esta escala. No utilizamos el código postal ya que no es un dato continuo.
subset_train <- na.omit(credit.Datos.Train)
subset_train$Deuda <- log10(subset_train$Deuda + 1)
subset_train$Anos_cotizados <- log10(subset_train$Anos_cotizados + 1)
subset_train$Calificacion_crediticia <- log10(subset_train$Calificacion_crediticia + 1)
subset_train$Ingresos <- log10(subset_train$Ingresos + 1)
pca <- subset_train[c("Edad", "Deuda", "Anos_cotizados", "Calificacion_crediticia", "Ingresos")]
Mediante un gráfico de dispersión vemos las distribuciones de las dos
primeras componentes calculadas diferenciando según
Aprobado.
pca_res <- prcomp(pca, scale = TRUE)
pca_df = as.data.frame(pca_res$x,stringsAsFactors=F)
pca_df = cbind(Aprobado = subset_train$Aprobado,pca_df)
ggplot(pca_df, aes(PC1, PC2)) +
modelr::geom_ref_line(h = 0) +
modelr::geom_ref_line(v = 0) +
geom_point(aes(color = Aprobado), size = 2, position='jitter',alpha = 0.4) +
xlab("Primera Componente Principal") +
ylab("Segunda Componente Principal") +
ggtitle("Dos primeras componentes principales")
Como vemos, la componente PC1 parece un buen predictor. Podemos visualizar su dispersión en el siguiente gráfico y utilizar su matriz de confusión para ver su rendimiento como predictor.
ggplot(pca_df, aes(x = PC1, fill = Aprobado)) +
geom_density(alpha = 0.5) +
labs(title = "Distribución de PC1 según Aprobado",
x = "PC1",
y = "Densidad") +
theme_minimal()
pca_df$Prediccion_simple <- ifelse(pca_df$PC1 > 0.2, "Aprobado", "Rechazado")
conf_mat <- table(Real = pca_df$Aprobado, Predicho = pca_df$Prediccion_simple)
conf_df <- as.data.frame(conf_mat)
conf_df <- conf_df %>%
mutate(Tipo = case_when(
Real == "Aprobado" & Predicho == "Aprobado" ~ "Verdadero Positivo",
Real == "Rechazado" & Predicho == "Aprobado" ~ "Falso Positivo",
Real == "Rechazado" & Predicho == "Rechazado" ~ "Verdadero Negativo",
Real == "Aprobado" & Predicho == "Rechazado" ~ "Falso Negativo"
))
ggplot(conf_df, aes(x = Predicho, y = Real, fill = Freq)) +
geom_tile(color = "white") +
geom_text(aes(label = paste0(Freq, "\n", Tipo)), size = 4) +
scale_fill_gradient(low = "white", high = "steelblue") +
labs(title = "Matriz de confusión para predictor basado en PC1",
x = "Predicción", y = "Clase real", fill = "Frecuencia") +
theme_minimal()
Este predictor acierta en el 75% de los casos.
De la componente PC1 podemos saber también cómo se calcula a partir de los campos numéricos originales:
loadings <- pca_res$rotation
loadings[,1]
## Edad Deuda Anos_cotizados
## -0.3497941 -0.4057191 -0.5065381
## Calificacion_crediticia Ingresos
## -0.5505842 -0.3915511
PC1 da más o menos la misma importancia a todas las variables,
relacionándose de forma negativa con todas. La que tiene un mayor
coeficiente es la Calificacion_crediticia. Cabe mencionar
que estos coeficientes se calculan sobre una distribución reescalada
para seguir una distribución normal, por lo que es correcto el análisis
de correlacionar una mayor importancia en la clasificación con respecto
al escalar dado por la tabla anterior.
Las cinco componentes PCX aparecen ordenadas según la proporción de la varianza explicada (PVE). Podemos verlo para estos datos.
VE <- pca_res$sdev^2
PVE <- VE / sum(VE)
round(PVE,2)
## [1] 0.35 0.22 0.17 0.15 0.11
La primera componente explica, pues, más de un tercio de la varianza. Junto con la segunda, consiguen explicar prácticamente el 60%.
Terminados el análisis y el preprocesado de los datos pasamos al entrenamiento de los modelos.
Seleccionamos la variable que vamos a usar como salida (Aprobado) y el resto se utilizarán como entrada. Vamos a utilizar la versión del dataset que hemos definido en el apartado de preprocesado.
credit_no_na <- na.omit(credit.Datos.Train)
credit.Var.Salida.Usada <- c("Aprobado")
credit.Vars.Entrada.Usadas <- setdiff(names(credit.Datos.Train), credit.Var.Salida.Usada)
Los modelos que utilicemos tienen que poder usar tanto de los atributos numéricos como los categóricos, aunque no necesariamente los dos tipos a la vez. Incluimos:
rpart: sencillo y basado en árboles de decisión
ranger: complejo y basado en bosques aleatorios
gbm: complejo y basado en descenso de gradiente
nnet: complejo y dentro del paradigma MLP
En dos de los modelos mencionados (GBM y MLP) hemos utilizado un grid
de hiperparámetros, a través del argumento tuneGrid para
encontrar las mejores combinaciones de valores para estos. Con el
comando expand.grid() se crea el data frame con las todas
las combinaciones a probar. Una vez se han entrenado los modelos,
podemos repetir el proceso con una estimacón mas afinada en base a los
resultados obtenidos. Finalmente nos quedamos con el modelo que obtiene
mejores resultados.
Como primera aproximación a la tarea de clasificación, nos pareció interesante usar una técnica sencilla como lo son los árboles de decisión. Este tipo de modelo predictivo recibe su nombre por tener una estructura de árbol: para clasificar un individuo se parte desde una raíz. En cada nodo, según el valor de una determinado atributo, se divide el conjunto de datos en dos o más ramas. Este proceso se repite hasta alcanzar una hoja, a la cual se asigna una única categoría. Las características por las que nos interesamos en este método son:
Las divisiones pueden realizarse en base a atributos tanto numéricos como categóricos, por lo que es capaz de aprovechar los patrones en variables categóricas sin necesidad de conversión.
La fácil interpretabilidad de los resultados obtenidos, que puede guiar nuestra elección de otros modelos a considerar.
Como punto negativo, esta técnica no es muy potente y es propensa al sobreajuste. Sin embargo, el objetivo de este ataque, más que obtener una solución final, es informar nuestra posterior elección de modelos, además de establecer unos resultados de clasificación base con los que contrastar los que obtendremos con otros modelos.
library(rpart)
set.seed(155)
tree_model <- train(
credit_no_na[credit.Vars.Entrada.Usadas],
credit_no_na[[credit.Var.Salida.Usada]],
method = "rpart",
metric = "Accuracy",
tuneLength = 10
)
tree_model
## CART
##
## 516 samples
## 16 predictor
## 2 classes: 'Aprobado', 'Rechazado'
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 516, 516, 516, 516, 516, 516, ...
## Resampling results across tuning parameters:
##
## cp Accuracy Kappa
## 0.00000000 0.8279908 0.6531977
## 0.07943262 0.8713235 0.7439150
## 0.15886525 0.8713235 0.7439150
## 0.23829787 0.8713235 0.7439150
## 0.31773050 0.8713235 0.7439150
## 0.39716312 0.8713235 0.7439150
## 0.47659574 0.8713235 0.7439150
## 0.55602837 0.8713235 0.7439150
## 0.63546099 0.8583286 0.7128778
## 0.71489362 0.6532674 0.2610764
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was cp = 0.5560284.
Hemos hecho uso de rpart con un parámetro
tuneLength de 10. Hemos probado distintos valores de este
hiperparámetro, pero los resultados que se obtienen son casi siempre los
mismos. Con este método, el único hiperpárametro a modificar es
cp (controla la complejidad del árbol, evitando sobreajuste
al requerir una mejora mínima en el error para permitir una nueva
división), que toma valores entre 0 y 1. Pues bien, para todo valor
entre 0.018 y 0.714, el árbol que se obtiene es exactamente el
mismo.
library(rpart.plot)
rpart.plot(tree_model$finalModel, type = 2, extra = 104, fallen.leaves = TRUE, main = "Árbol de Decisión")
Evaluando los resultados obtenidos, nos damos cuenta que el árbol encontrado con mejor accuracy no es otro que el decisor que ya habíamos considerado en la fase de análisis exploratorio: si existe impago previo se rechaza la solicitud, si no, se aprueba.
El siguiente modelo que vamos a utilizar es Random Forest. Los Bosques Aleatorios son una técnica de aprendizaje automático que opera construyendo múltiples árboles de decisión durante el entrenamiento y generando la clase que es elegida con mayor frecuencia en todos los árboles individuales para la clasificación.
En concreto usamos ranger que está dentro de la librería caret. La razón principal para elegir ranger es su eficiencia. ranger es conocida por ser una de las implementaciones más rápidas de Bosques Aleatorios disponibles en R, o que nos permitirá entrenar modelos de Bosques Aleatorios en grandes conjuntos de datos de manera eficiente.
Los dos hiperparámetros más influyentes para optimizar el rendimiento del modelo son mtry y splitrule:
mtry: Controla el número de variables predictoras que se consideran aleatoriamente en cada división de un nodo al construir un árbol individual.
splitrule: Define el criterio utilizado para decidir cómo se realiza la mejor división en cada nodo de un árbol. Las opciones disponibles para esta parámetro son:
gini (por defecto): Utiliza el índice de Gini para medir la impureza.
extratrees: Implementa el algoritmo de “Extremely Randomized Trees”. En lugar de buscar la mejor división para cada variable (como lo hace Gini), elige puntos de división aleatorios para cada variable y luego selecciona la mejor de estas divisiones aleatorias.
Además, utilizando un tuneLength de 15 ya podemos comprobar todas las combinaciones posibles de hiperparámetros en nuestro caso. Para el parámetro num.trees hemos probado con valores desde 50 hasta 1000 y no hemos obtenido diferencias notables en los valores en este rango, por lo que decidimos utilizar 50 ya que el tiempo de computación es mucho menor.
Realizamos el entrenamiento:
n_reps=3
n_folds=10
set.seed(68)
fitControl <- trainControl(
method = "repeatedcv",
number = n_folds,
repeats = n_reps,
verboseIter = FALSE,
)
credit.ranger <- train(
credit_no_na[credit.Vars.Entrada.Usadas],
credit_no_na[[credit.Var.Salida.Usada]],
method = "ranger",
verbose = FALSE,
tuneLength = 15,
trControl = fitControl,
num.trees = 50,
importance = 'impurity'
)
credit.ranger
## Random Forest
##
## 516 samples
## 16 predictor
## 2 classes: 'Aprobado', 'Rechazado'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 465, 465, 464, 465, 464, 465, ...
## Resampling results across tuning parameters:
##
## mtry splitrule Accuracy Kappa
## 2 gini 0.8728929 0.7427871
## 2 extratrees 0.8773307 0.7524617
## 3 gini 0.8799456 0.7579299
## 3 extratrees 0.8857772 0.7698214
## 4 gini 0.8811899 0.7609014
## 4 extratrees 0.8740763 0.7469396
## 5 gini 0.8779597 0.7545235
## 5 extratrees 0.8780337 0.7547231
## 6 gini 0.8811406 0.7612200
## 6 extratrees 0.8773559 0.7534241
## 7 gini 0.8785760 0.7560194
## 7 extratrees 0.8799079 0.7585862
## 8 gini 0.8740632 0.7469176
## 8 extratrees 0.8824473 0.7636032
## 9 gini 0.8792431 0.7576588
## 9 extratrees 0.8799456 0.7586401
## 10 gini 0.8799200 0.7587888
## 10 extratrees 0.8793162 0.7571140
## 11 gini 0.8786254 0.7563495
## 11 extratrees 0.8779099 0.7544731
## 12 gini 0.8774676 0.7536245
## 12 extratrees 0.8806239 0.7602259
## 13 gini 0.8780216 0.7547268
## 13 extratrees 0.8747285 0.7482210
## 14 gini 0.8792920 0.7574888
## 14 extratrees 0.8857148 0.7703562
## 15 gini 0.8773680 0.7535782
## 15 extratrees 0.8805736 0.7598006
## 16 gini 0.8793176 0.7577196
## 16 extratrees 0.8824971 0.7636341
##
## Tuning parameter 'min.node.size' was held constant at a value of 1
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were mtry = 3, splitrule = extratrees
## and min.node.size = 1.
varImp(credit.ranger)
## ranger variable importance
##
## Overall
## Impago_previo 100.0000
## Calificacion_crediticia 23.1176
## Esfuerzo_financiero 21.4281
## Trabaja 19.7108
## Anos_cotizados 18.9877
## Codigo_postal 12.4039
## Ingresos 12.0913
## Deuda 10.9468
## Nivel_educativo 8.8516
## Edad 8.7345
## Etnia 5.7209
## Licencia_de_conducir 2.6469
## Sexo 1.2724
## Ciudadano 0.2462
## Estado_civil 0.1240
## Es_cliente 0.0000
Obtenemos una accuracy de 88.57, utilizando mtry=3 y splitrule = extratrees.
También nos damos cuenta que, como ocurría en rpart, la variable de mayor importancia es Impago previo. También le da bastante importancia a la variable añadida por nosotros (Esfuerzo_financiero), lo cual nos hace pensar que ha sido una buena inclusión.
Gradient boosting (o potenciación del gradiente en castellano) funciona construyendo un conjunto de modelos predictivos que ejecutan decisiones simples sobre los datos. Estos decisores son normalmente árboles de decisión, aunque también puede emplear métodos de regresión lineal, splines o redes neuronales poco profundas, entre otros.
La justificación en el uso de árboles radica en su flexibilidad tratando todo tipo de datos: numéricos y categóricos, capacidad para capturar interacciones entre variables, eficiencia computacional e interpretabilidad.
Gbm produce nuevos modelos de forma secuencial, buscando minimizar en cada iteración una función de pérdida arbitraria que mide la diferencia entre las predicciones y los valores reales. El nuevo modelo se ajusta a los “pseudo-residuales”, que son los gradientes negativos de la función de pérdida con respecto a las predicciones actuales. Este proceder es lo que le da el nombre de “gradient”, pues lo que se hace en esencia es minimizar en la dirección que más decrece el gradiente de la función de pérdidas.
Escogemos gbm ya que:
Permite trabajar con combinaciones de datos numéricos y categóricos.
Captación de relaciones entre componentes, tal cual lo hacen otras técnicas que usan árboles de decisión.
Perfecto para clasificación binaria con la función de pérdida de desviación binaria (por defecto).
En primer lugar, vamos a definir un grid de hiperparámetros para encontrar las condiciones ideales con las que utilizar gbm. Vamos a hacer una exploración para el número de árboles (n.trees), para la profundidad máxima de los mismos (interaction.depth), la tasa de aprendizaje (shrinkage) y mínimo de observaciones por nodo (n.minobsinnode). También fijamos una semilla para poder recrear los resultados.
Para la creación del grid seguimos el método de ajustar el grano. En primer lugar empezamos con valores de n.trees en un mayor rango (100, 200, 300). Nos dimos cuenta que funcionaba mejor con valores más bajos, por lo que ajustamos a (25,50,75). Como el que mejor estaba funcionando era 50, ajustamos a (45, 50, 55) que fueron los valores que finalmente utilizamos. De forma similar para los otros tres parámetros.
library(caret)
set.seed(1234)
gbm_grid <- expand.grid(
n.trees = c(45, 50, 55),
interaction.depth = c(3,6,7),
shrinkage = c(0.09, 0.1, 0.12),
n.minobsinnode = c(10, 15,20)
)
n_folds <- 10
n_reps <- 3
Posteriormente, definimos el control de entrenamiento para determinar cómo realizar la validación cruzada y qué metricas usar para evaluar la combinación de los hiperparámetros. Vamos a utilizar validación cruzada repetitiva con 10 folds y 3 repeticiones.
fitControl <- trainControl(
method = "repeatedcv",
number = n_folds,
repeats = n_reps,
allowParallel = TRUE
)
Finalmente, realizamos el entrenamiento del modelo, especificando la columna de salida “Aprobado”.
gbm_modelo <- train(
credit_no_na[credit.Vars.Entrada.Usadas],
credit_no_na[[credit.Var.Salida.Usada]],
method = "gbm",
trControl = fitControl,
tuneGrid = gbm_grid,
verbose = FALSE,
bag.fraction = 0.75
)
gbm_modelo
## Stochastic Gradient Boosting
##
## 516 samples
## 16 predictor
## 2 classes: 'Aprobado', 'Rechazado'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 464, 465, 464, 464, 464, 465, ...
## Resampling results across tuning parameters:
##
## shrinkage interaction.depth n.minobsinnode n.trees Accuracy Kappa
## 0.09 3 10 45 0.8714145 0.7414575
## 0.09 3 10 50 0.8713898 0.7411726
## 0.09 3 10 55 0.8707362 0.7398796
## 0.09 3 15 45 0.8687996 0.7358640
## 0.09 3 15 50 0.8726337 0.7433486
## 0.09 3 15 55 0.8707353 0.7394551
## 0.09 3 20 45 0.8668645 0.7319577
## 0.09 3 20 50 0.8688001 0.7357091
## 0.09 3 20 55 0.8681339 0.7342962
## 0.09 6 10 45 0.8719554 0.7420917
## 0.09 6 10 50 0.8693908 0.7370076
## 0.09 6 10 55 0.8706855 0.7393270
## 0.09 6 15 45 0.8764929 0.7510608
## 0.09 6 15 50 0.8784411 0.7548993
## 0.09 6 15 55 0.8777875 0.7535700
## 0.09 6 20 45 0.8739530 0.7460250
## 0.09 6 20 50 0.8726458 0.7433068
## 0.09 6 20 55 0.8733245 0.7447973
## 0.09 7 10 45 0.8784662 0.7550671
## 0.09 7 10 50 0.8797106 0.7574232
## 0.09 7 10 55 0.8804019 0.7587960
## 0.09 7 15 45 0.8706985 0.7391930
## 0.09 7 15 50 0.8681218 0.7341547
## 0.09 7 15 55 0.8726342 0.7431256
## 0.09 7 20 45 0.8784034 0.7551161
## 0.09 7 20 50 0.8765171 0.7511697
## 0.09 7 20 55 0.8771707 0.7523574
## 0.10 3 10 45 0.8687880 0.7359302
## 0.10 3 10 50 0.8694165 0.7371584
## 0.10 3 10 55 0.8687880 0.7355993
## 0.10 3 15 45 0.8694406 0.7371854
## 0.10 3 15 50 0.8713763 0.7411704
## 0.10 3 15 55 0.8707353 0.7398139
## 0.10 3 20 45 0.8713395 0.7408489
## 0.10 3 20 50 0.8706855 0.7393703
## 0.10 3 20 55 0.8738906 0.7458279
## 0.10 6 10 45 0.8745945 0.7476502
## 0.10 6 10 50 0.8726714 0.7436204
## 0.10 6 10 55 0.8752481 0.7487790
## 0.10 6 15 45 0.8751978 0.7485347
## 0.10 6 15 50 0.8758388 0.7494359
## 0.10 6 15 55 0.8771590 0.7523927
## 0.10 6 20 45 0.8765804 0.7512150
## 0.10 6 20 50 0.8772214 0.7525768
## 0.10 6 20 55 0.8759017 0.7496463
## 0.10 7 10 45 0.8752732 0.7487906
## 0.10 7 10 50 0.8752983 0.7487677
## 0.10 7 10 55 0.8740037 0.7459925
## 0.10 7 15 45 0.8739283 0.7455662
## 0.10 7 15 50 0.8720173 0.7418818
## 0.10 7 15 55 0.8687745 0.7353436
## 0.10 7 20 45 0.8758514 0.7499555
## 0.10 7 20 50 0.8764673 0.7510515
## 0.10 7 20 55 0.8784155 0.7548864
## 0.12 3 10 45 0.8714019 0.7412620
## 0.12 3 10 50 0.8733250 0.7450175
## 0.12 3 10 55 0.8719810 0.7420966
## 0.12 3 15 45 0.8700444 0.7381897
## 0.12 3 15 50 0.8732998 0.7446218
## 0.12 3 15 55 0.8694034 0.7368385
## 0.12 3 20 45 0.8687754 0.7352955
## 0.12 3 20 50 0.8707106 0.7392561
## 0.12 3 20 55 0.8726467 0.7431383
## 0.12 6 10 45 0.8836819 0.7656035
## 0.12 6 10 50 0.8830163 0.7642309
## 0.12 6 10 55 0.8804517 0.7589838
## 0.12 6 15 45 0.8713144 0.7403951
## 0.12 6 15 50 0.8713647 0.7402485
## 0.12 6 15 55 0.8668765 0.7313605
## 0.12 6 20 45 0.8784909 0.7549484
## 0.12 6 20 50 0.8764803 0.7509522
## 0.12 6 20 55 0.8739153 0.7457079
## 0.12 7 10 45 0.8739036 0.7458889
## 0.12 7 10 50 0.8739418 0.7459762
## 0.12 7 10 55 0.8758770 0.7498115
## 0.12 7 15 45 0.8746322 0.7472115
## 0.12 7 15 50 0.8758765 0.7493578
## 0.12 7 15 55 0.8752602 0.7482628
## 0.12 7 20 45 0.8764924 0.7512178
## 0.12 7 20 50 0.8771460 0.7523710
## 0.12 7 20 55 0.8765301 0.7510315
##
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were n.trees = 45, interaction.depth =
## 6, shrinkage = 0.12 and n.minobsinnode = 10.
El resultado óptimo de accuracy (0.8836819) se ha obtenido con los siguiente parámetros:
n.trees = 45
interaction.depth = 6
shrinkage = 0.12
n.minobsinnode = 10
Al dar como resultado una profundidad de 6, sabemos que el árbol final de decisión utiliza seis comprobaciones para la clasificación. Podemos ver la evolución de la acurracy en función de la tasa de aprendizaje y el resto de variables con la siguientes gráficas:
plot(gbm_modelo, metric = "Accuracy", plotType = "scatter", ylim = c(0.85, 0.89))
En lo relativo a la importancia de las variables, podemos observarla ejecutando el siguiente comando:
varImp(gbm_modelo)
## gbm variable importance
##
## Overall
## Impago_previo 100.0000
## Nivel_educativo 17.9465
## Codigo_postal 11.2446
## Esfuerzo_financiero 10.8537
## Deuda 6.6734
## Ingresos 5.7575
## Edad 3.6635
## Trabaja 3.0078
## Anos_cotizados 1.8489
## Calificacion_crediticia 1.6783
## Etnia 1.5437
## Ciudadano 0.6836
## Sexo 0.6345
## Estado_civil 0.6119
## Licencia_de_conducir 0.1393
## Es_cliente 0.0000
Nuevamente, el impago previo es la variable más determinante. Sorprendentemente, Nivel_educativo es el segundo con mayor importancia. También cabe destacar que Esfuerzo_financiero, la variable que hemos creado nosotros es una de las que tiene mayor relevancia.
Finalmente, queremos entrenar también un perceptrón multicapa (MLP)
para predecir la variable Aprobado. Para esto hacemos uso
de la librería nnet.
Un MLP es una red neuronal compuesta por varias capas de neuronas conectadas entre sí:
Una capa de entrada, por la que se introducen los datos al modelo. Cada neurona de entrada representa una variable o característica del conjunto de datos.
Una o varias capas ocultas, que son las responsables del aprendizaje interno del modelo. Las neuronas de cada capa oculta están conectadas a todas las neuronas de la capa anterior.
Una capa de salida, que genera la predicción final. En este caso, puesto que la nuestra es una tarea de clasificación binaria, esta esta formada por una única neurona con función sigmoide.
Dado que un MLP sólo procesa datos númericos y debido a la gran
importancia del campo Impago_previo, decidimos añadir un
nuevo atributo al dataset, llamado Impago_previo_num. Este
campo es sencillamente una transformación de Impago_previo
a un campo numérico, para poder pasárselo como entrada a nuestro
MLP.
credit_copy <- credit
credit_copy$Esfuerzo_financiero <- credit_copy$Deuda + credit_copy$Anos_cotizados + abs(credit_copy$Ingresos-1)*0.375
credit_copy$Impago_previo_num <- ifelse(credit_copy$Impago_previo == "f", 1, 0)
campos_numericos_nuevo <- c(campos_numericos, 18)
campos_numericos_nuevo
## [1] 2 3 8 11 14 15 17 18
# Separación de conjuntos
credit_copy_train <- credit_copy[credit.trainIdx, ]
credit_copy_test <- credit_copy[-credit.trainIdx, ]
A continuación, elegimos los campos que usaremos como entrada para el
modelo. Probamos un primer preprocesado eliminando el
Codigo_postal de las variables a introducir al MLP. Sin
embargo nos sorprendió que se producía una ligera disminución en
Accuracy. Esperábamos que este no aportara ningún valor
como campo numérico. No obstante, dado que parece producir una pequeña
mejora, lo hemos conservado, siendo este el segundo preprocesado.
mlp_inputs <- campos_numericos_nuevo
campos_numericos_nuevo
## [1] 2 3 8 11 14 15 17 18
mlp_inputs_no_cp <- c(2, 3, 8, 11, 15, 17, 18)
El ajuste de hiperparámetros lo llevaremos a cabo usando
tuneGrid, y modificaremos los valores de size
y decay:
size hace referencia al número de neuronas en la
capa oculta de la red. El conjunto de valores que hemos elegido es
[3, 7, 11, 15]
Weight decay es una técnica para evitar el sobrentrenamiento.
Consiste en minimizar los pesos de las neuronas junto con el error de
clasificación cometido. decay controla la intensidad de la
regularización. Los valores elegidos han sido
[0, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5, 1]
credit_no_na <- credit_copy_train
#tomamos logaritmo de variables en campos_log
credit_log <- credit_no_na
credit_log[campos_log] <- log(credit_log[campos_log] + 1)
#normalizamos
campos_numericos
## [1] 2 3 8 11 14 15 17
credit_scaled <- credit_log
summary(credit_scaled)
## Sexo Edad Deuda Estado_civil Es_cliente
## a :167 Min. :13.75 Min. : 0.000 l : 2 g :419
## b :374 1st Qu.:22.50 1st Qu.: 1.000 u :419 gg : 2
## NA's: 12 Median :27.83 Median : 2.750 y :126 p :126
## Mean :31.32 Mean : 4.795 t : 0 NA's: 6
## 3rd Qu.:37.33 3rd Qu.: 7.500 NA's: 6
## Max. :80.25 Max. :28.000
## NA's :12
## Nivel_educativo Etnia Anos_cotizados Impago_previo Trabaja
## c :103 v :321 Min. :0.0000 f:269 f:319
## q : 58 h :101 1st Qu.:0.1527 t:284 t:234
## i : 53 bb : 48 Median :0.6931
## w : 48 ff : 48 Mean :0.8114
## ff : 46 j : 7 3rd Qu.:1.2528
## (Other):236 (Other): 19 Max. :3.3844
## NA's : 9 NA's : 9
## Calificacion_crediticia Licencia_de_conducir Ciudadano Codigo_postal
## Min. :0.0000 f:299 g:502 Min. : 0.0
## 1st Qu.:0.0000 t:254 p: 8 1st Qu.: 73.0
## Median :0.0000 s: 43 Median : 160.0
## Mean :0.6891 Mean : 183.5
## 3rd Qu.:1.3863 3rd Qu.: 280.0
## Max. :4.2195 Max. :2000.0
## NA's :13
## Ingresos Aprobado Esfuerzo_financiero Impago_previo_num
## Min. : 0.000 Aprobado :246 Min. : 0.000 Min. :0.0000
## 1st Qu.: 0.000 Rechazado:307 1st Qu.: 4.625 1st Qu.:0.0000
## Median : 1.792 Median : 13.875 Median :0.0000
## Mean : 3.019 Mean : 413.702 Mean :0.4864
## 3rd Qu.: 5.994 3rd Qu.: 160.125 3rd Qu.:1.0000
## Max. :11.513 Max. :37521.625 Max. :1.0000
##
credit_scaled[campos_numericos_nuevo] <- scale(credit_scaled[campos_numericos_nuevo])
sum(is.na(credit_scaled))
## [1] 67
set.seed(123)
ctrl <- trainControl(method = "cv", number = 10)
mlp_cv <- train(
credit_scaled[mlp_inputs],
credit_scaled[["Aprobado"]],
method = "nnet",
trControl = ctrl,
tuneGrid = expand.grid(
size = c(3, 7, 11, 15),
decay = c(0, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5, 1)
),
maxit = 200,
trace = FALSE
)
print(mlp_cv)
## Neural Network
##
## 553 samples
## 8 predictor
## 2 classes: 'Aprobado', 'Rechazado'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 498, 497, 498, 499, 497, 498, ...
## Resampling results across tuning parameters:
##
## size decay Accuracy Kappa
## 3 0.000 0.8516138 0.7002978
## 3 0.003 0.8641974 0.7257857
## 3 0.010 0.8565631 0.7105817
## 3 0.030 0.8750984 0.7484643
## 3 0.100 0.8713003 0.7410661
## 3 0.300 0.8773441 0.7534595
## 3 0.500 0.8773441 0.7536109
## 3 1.000 0.8716461 0.7422107
## 7 0.000 0.8487889 0.6934416
## 7 0.003 0.8386575 0.6742771
## 7 0.010 0.8509248 0.6981161
## 7 0.030 0.8636750 0.7254501
## 7 0.100 0.8769634 0.7518704
## 7 0.300 0.8810022 0.7607830
## 7 0.500 0.8755387 0.7501206
## 7 1.000 0.8734602 0.7462522
## 11 0.000 0.8461253 0.6901444
## 11 0.003 0.8259090 0.6472568
## 11 0.010 0.8426883 0.6805930
## 11 0.030 0.8548255 0.7064143
## 11 0.100 0.8540653 0.7056813
## 11 0.300 0.8789041 0.7566262
## 11 0.500 0.8813080 0.7619020
## 11 1.000 0.8755010 0.7503457
## 15 0.000 0.8040528 0.6040669
## 15 0.003 0.8172471 0.6315260
## 15 0.010 0.8047214 0.6040123
## 15 0.030 0.8445849 0.6859169
## 15 0.100 0.8675734 0.7331494
## 15 0.300 0.8847864 0.7684581
## 15 0.500 0.8755387 0.7501206
## 15 1.000 0.8755010 0.7503457
##
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were size = 15 and decay = 0.3.
plot(mlp_cv)
Observando los resultados obtenidos podemos comprobar que el efecto
de decay es notable y muy positivo. Al ser un conjunto de
datos pequeño, existe una tendencia fuerte a sobreajustar. Forzar los
pesos a tomar valores pequeños limita este efecto.
También vemos que los resultados obtenidos son similares independientemente del número de neuronas en la capa oculta. De este hecho deducimos que este dataset no contiene patrones que se puedan explotar fácilmente para mejorar la clasificación, más allá de alguna tenue mejora sobre los patrones superficiales ya encontrados por nosotros.
Los mejores resultados se obtienen con un valor de size
15, con un decay de 0.3. En realidad este máximo parece
causa de ligeras desviaciones aleatorias más que de una adecuación mayor
del size elegido, lo que se puede deducir del hecho de que
las cuatro gráficas se entrecruzan constantemente.
En esta sección llevamos a cabo la comparación de los cuatro modelos.
Todos los modelos se han entrenado bajo una misma configuración de validación cruzada (10 particiones repetidas 3 veces), utilizando como métrica de evaluación la exactitud (Accuracy). El objetivo es identificar qué modelo se ajusta mejor a los datos y generaliza con mayor eficacia.
set.seed(1234)
credit.trainCtrl.3cv10.resampAll <- trainControl(
method = "repeatedcv"
,number = 10 ,repeats = 3,
verboseIter=F,
returnResamp = "all"
)
credit.comparison.rpart<- suppressWarnings(train(credit.Datos.Train[credit.Vars.Entrada.Usadas],
credit.Datos.Train[[credit.Var.Salida.Usada]],
method='rpart',
metric = "Accuracy",
tuneLength = 10,
trControl = credit.trainCtrl.3cv10.resampAll
))
credit.comparison.gbm <- suppressWarnings(train(
credit_no_na[credit.Vars.Entrada.Usadas],
credit_no_na[[credit.Var.Salida.Usada]],
method = 'gbm',
tuneGrid = gbm_grid,
verbose = FALSE,
trControl = credit.trainCtrl.3cv10.resampAll,
bag.fraction = 0.75
))
credit.comparison.ranger <- suppressWarnings(train(
credit_no_na[credit.Vars.Entrada.Usadas],
credit_no_na[[credit.Var.Salida.Usada]],
method = 'ranger',
verbose = FALSE,
metric = "Accuracy",
tuneLength = 15,
trControl = credit.trainCtrl.3cv10.resampAll
))
credit.comparison.mlp <- suppressWarnings(train(
credit_scaled[mlp_inputs],
credit_scaled[["Aprobado"]],
method = "nnet",
trControl = credit.trainCtrl.3cv10.resampAll,
tuneGrid = expand.grid(
size = c(3, 5, 7, 9, 11),
decay = c(0, 0.003, 0.01, 0.03, 0.1, 0.3, 1)
),
maxit = 200,
trace = FALSE
))
models <- list(rpart = credit.comparison.rpart, gbm = gbm_modelo , ranger = credit.comparison.ranger, mlp = credit.comparison.mlp)
credit.rsamp.comparison <- resamples(models)
## Warning in resamples.default(models): 'rpart' did not have
## 'returnResamp="final"; the optimal tuning parameters are used
## Warning in resamples.default(models): 'ranger' did not have
## 'returnResamp="final"; the optimal tuning parameters are used
## Warning in resamples.default(models): 'mlp' did not have 'returnResamp="final";
## the optimal tuning parameters are used
summary(credit.rsamp.comparison)
##
## Call:
## summary.resamples(object = credit.rsamp.comparison)
##
## Models: rpart, gbm, ranger, mlp
## Number of resamples: 30
##
## Accuracy
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## rpart 0.7962963 0.8363636 0.8624579 0.8646946 0.8909091 0.9444444 0
## gbm 0.7884615 0.8503017 0.8846154 0.8836819 0.9215686 0.9615385 0
## ranger 0.7678571 0.8424272 0.8808081 0.8799655 0.9103084 0.9636364 0
## mlp 0.7843137 0.8518519 0.8774929 0.8737968 0.8991979 1.0000000 0
##
## Kappa
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## rpart 0.6024096 0.6760434 0.7265795 0.7311502 0.7827515 0.8888889 0
## gbm 0.5705706 0.6966281 0.7692308 0.7656035 0.8420909 0.9230769 0
## ranger 0.5210526 0.6857595 0.7584301 0.7563011 0.8190519 0.9267643 0
## mlp 0.5661253 0.7049180 0.7555296 0.7461759 0.7943952 1.0000000 0
dotplot(credit.rsamp.comparison, metric = "Accuracy",
scales =list(x = list(relation = "free")))
Tras analizar los resultados obtenidos por cada modelo en validación cruzada, se observa que tanto el modelo gbm (Gradient Boosting Machine) como el ranger obtienen resultados muy parecidos, siendo los de gbm ligeramente superiores, tanto en el primer cuartil como en la media. En comparación, el modelo rpart, aunque interpretable, ha mostrado un rendimiento inferior, como es habitual en árboles individuales. El modelo nnet (red neuronal) ha obtenido resultados competitivos, pero no han superado a los otros dos.
Por tanto, se selecciona el modelo gbm como el más adecuado para tratar los casos de test. Veamos los resultados que obtiene:
preds<-predict( gbm_modelo,
newdata=credit.Datos.Test[credit.Vars.Entrada.Usadas])
confusionMatrix(preds, credit.Datos.Test[[credit.Var.Salida.Usada]])
## Confusion Matrix and Statistics
##
## Reference
## Prediction Aprobado Rechazado
## Aprobado 54 10
## Rechazado 7 66
##
## Accuracy : 0.8759
## 95% CI : (0.8088, 0.926)
## No Information Rate : 0.5547
## P-Value [Acc > NIR] : 5.29e-16
##
## Kappa : 0.75
##
## Mcnemar's Test P-Value : 0.6276
##
## Sensitivity : 0.8852
## Specificity : 0.8684
## Pos Pred Value : 0.8438
## Neg Pred Value : 0.9041
## Prevalence : 0.4453
## Detection Rate : 0.3942
## Detection Prevalence : 0.4672
## Balanced Accuracy : 0.8768
##
## 'Positive' Class : Aprobado
##
Obtenemos una Accuracy del 87.5%, un resultado más que
aceptable, aunque nos damos cuenta que es inferior a lo que habíamos
obtenido en la validación con el subconjunto de Train.
En esta práctica hemos realizado un análisis de un dataset muy completo desde el punto de vista académico, ya que mantiene un equilibrio entre los tipos de atributos y requiere una reflexión sobre el tratamiento de outliers y valores nulos.
El análisis exploratorio de datos mostró una gran influencia del
campo Impago_previo en la clasificación de los individuos.
Por ello, tanto los modelos predictivos hechos a mano como los modelos
de aprendizaje automático han dado la mayor importancia a este
campo.
A pesar de que el mejor modelo de los utilizados es gbm, tanto en rendimiento medio como en consistencia, todos los modelos seleccionados tienen un rendimiento parecido, al menos en los datos de entrenamiento.
Como valoración general, estamos bastante contentos con el trabajo realizado y nos ha parecido una buena forma de poner en práctica los conocimientos de la asignatura. A pesar de estar sesgados por conocer qué representan los atributos, creemos que el tratamiento final de los datos ha sido riguroso. Con respecto a la utilización de los modelos, gracias a la exploración de meta-parámetros con tuneGrid hemos podido entender la influencia de estos, teniendo que comprobar primero la teoría para entenderlos.